Analiza długości śledzi

Długość śledzia

p_length <- ggplot(dataset, aes(length)) +
     geom_histogram(binwidth = 0.5) + 
     scale_x_discrete(name="Długość [cm]", limits= seq(min(dataset$length), max(dataset$length), by=1) ) +
     ylab("Liczba obserwacji") +
     ggtitle("Długość śledzia") + 
     theme(plot.title = element_text(hjust = 0.5))

Rozkład wartości długości śledzia jest zbliżony do rozkładu normalnego o średniej ok. 25.5 cm.

Zagęszczenie planktonu

p_cfin1 <- ggplot(dataset, aes(cfin1)) + geom_histogram(binwidth = 1.0) + xlab("Zagęszczenie planktonu") +
  ylab("Liczba obserwacji") + ggtitle("Calanus finmarchicus gat. 1") + theme(plot.title = element_text(hjust = 0.5))

p_cfin2 <- ggplot(dataset, aes(cfin2)) + geom_histogram(binwidth = 1.0) + xlab("Zagęszczenie planktonu") +  
  ylab("Liczba obserwacji") + ggtitle("Calanus finmarchicus gat. 2") + theme(plot.title = element_text(hjust = 0.5))

grid.arrange(p_cfin1, p_cfin2, ncol = 2, nrow = 1)

p_chel1 <- ggplot(dataset, aes(chel1)) +  geom_histogram(binwidth = 1.0) + xlab("Zagęszczenie planktonu") +
  ylab("Liczba obserwacji") + ggtitle("Calanus helgolandicus gat. 1") + theme(plot.title = element_text(hjust = 0.5))

p_chel2 <- ggplot(dataset, aes(chel2)) +  geom_histogram(binwidth = 1.0) + xlab("Zagęszczenie planktonu") +  
  ylab("Liczba obserwacji") + ggtitle("Calanus helgolandicus gat. 2") + theme(plot.title = element_text(hjust = 0.5))

grid.arrange(p_chel1, p_chel2, ncol = 2, nrow = 1)

p_lcop1 <- ggplot(dataset, aes(x=lcop1)) +  geom_histogram(binwidth = 1.0) + xlab("Zagęszczenie planktonu") +
  ylab("Liczba obserwacji") + ggtitle("Widłonogi gat. 1") + theme(plot.title = element_text(hjust = 0.5))

p_lcop2 <- ggplot(dataset, aes(x=lcop2)) +  geom_histogram(binwidth = 1.0) +   xlab("Zagęszczenie planktonu") +  
  ylab("Liczba obserwacji") + ggtitle("Widłonogi gat. 2") + theme(plot.title = element_text(hjust = 0.5))

grid.arrange(p_lcop1, p_lcop2, ncol = 2, nrow = 1)

Z wykresów wynika, że w większości przypadkóW zagęszczenie planktonu ma rozkład prawostronnie skośny. Dla Calanus finmarchicus gat. 1 oraz Widłonogów gat. 1 obserwujemy występowanie wartości odstających tzw. outlierów. Mogą być one skutkiem błędów popełnionych przy zbieraniu danych i negatywnie wpływać na dalszą analizę. Aby nie tracić zupełnie danych z wierszy, dla których niektóre atrybuty są outlierami, zastąpiono je wartościami NA, a następnie przeprowadzono interpolację.

replace_cfin1 <- which(dataset$cfin1 > 20)
replace_lcop1 <- which(dataset$cfin1 > 90)

for (i in replace_cfin1){
  dataset[i, "cfin1"] <- NA
}
for (j in replace_lcop1){
  dataset[i, "lcop1"] <- NA
}

dataset$cfin1 <- na.interpolation(dataset$cfin1)
dataset$lcop1 <- na.interpolation(dataset$lcop1)

Po przetworzeniu wartości odstających rozkład Calanus finmarchicus gat. 1 i Widłonogów gat. 1 jest bardziej równomierny i prezentuje się następująco:

Dane dotyczące liczby śledzi

p_fbar <- ggplot(dataset, aes(fbar)) + geom_histogram(binwidth = 0.05) + xlab("Ułamek pozostawionego narybku")+
  ylab("Liczba obserwacji") + ggtitle("Natężenie połowów w regionie") + theme(plot.title = element_text(hjust = 0.5))

p_recr <- ggplot(dataset, aes(recr)) + geom_histogram(binwidth = 50000.0) + xlab("Liczba śledzi")+
  ylab("Liczba obserwacji") + ggtitle("Roczny narybek") + theme(plot.title = element_text(hjust = 0.5))

p_cumf <- ggplot(dataset, aes(cumf)) + geom_histogram(binwidth = 0.02) + xlab("Ułamek pozostawionego narybku")+
  ylab("Liczba obserwacji") + ggtitle("Łączne roczne natężenie połowów") + theme(plot.title = element_text(hjust = 0.5))

p_totaln <- ggplot(dataset, aes(totaln)) + geom_histogram(binwidth = 50000.0) + xlab("Liczba śledzi")+
  ylab("Liczba obserwacji") + ggtitle("Łączna liczba złowionych ryb") + theme(plot.title = element_text(hjust = 0.5))

grid.arrange(p_fbar, p_recr, p_cumf, p_totaln, ncol = 2, nrow = 2)

Dane dotyczące liczby śledzi w łowiskach mają rozkład dużo bardziej równomierny niż dane dotyczące zagęszczenia planktonu. Natężenia połowów mają rozkład wielomodalny.

Warunki przyrodnicze

p_sst <- ggplot(dataset, aes(sst)) + geom_histogram(binwidth = 0.1) + xlab("Temperatura [°C]")+ 
  ylab("Liczba obserwacji")+ ggtitle("Temperatura przy powierzchni wody")+ theme(plot.title = element_text(hjust=0.5))

p_sal <- ggplot(dataset, aes(sal)) + geom_histogram(binwidth = 0.01) + xlab("Poziom zasolenia [Knudsen ppt]")+ 
  ylab("Liczba obserwacji") + ggtitle("Zasolenie wody") + theme(plot.title=element_text(hjust=0.5))

p_xmonth <- ggplot(dataset, aes(xmonth)) + geom_histogram(binwidth = 1.0) + xlab("Numer miesiąca")+
  ylab("Liczba obserwacji") + ggtitle("Miesiąc połowu") +  theme(plot.title=element_text(hjust=0.5))

p_nao <- ggplot(dataset, aes(nao)) + geom_histogram(binwidth = 0.2) + xlab("Oscylacja [mb]")+ 
  ylab("Liczba obserwacji") + ggtitle("Oscylacja północnoatlantycka") + theme(plot.title=element_text(hjust=0.5))

grid.arrange(p_sst, p_sal, p_xmonth, p_nao, ncol = 2, nrow = 2)

Najczęściej odnotowaną temperaturą wody było 13.6°C. Może być to związane z faktem, że najwięcej danych zostało zebranych w miesiącach letnich - szczególnie w sierpniu. Zasolenie wody koncentruje się głównie w przedziale 35.5 - 35.55 z niewielkimi odstępstwami.

Przygotowanie danych

Aby skonstruować dobry model uczenia maszynowego, konieczne jest wcześniejsze uważne przygotowanie danych. Wśród atrybutów znajduje się kolumna xmonth oznaczająca numer miesiąca. Mimo tego, że wartości reprezentowane są jako liczby, jest to atrybut kategoryczny. Dlatego jego wartości zostały zamienione za pomocą one-hot encoding - dodano 12 kolumn, z których dla każdego wiersza tylko jedna przyjmuje wartość 1, a pozostałe - 0.

cols <- colnames(dataset)
dataset2 <- cbind(dataset, dummy(dataset$xmonth, sep = "_"))
cols <- c(cols, c("jan", "feb", "mar", "apr", "may", "jun", "jul", "aug", "sep", "oct", "nov", "dec"))
colnames(dataset2) <- cols
dataset2['xmonth'] <- NULL
head(dataset2) %>%
  kable() %>%
  kable_styling() %>%
  scroll_box(width = "900px")

Kolejnym krokiem jest podział zbioru na zbiór uczący, testowy i walidacyjny. Przyjęto proporcje odpowiednio: 7:2:1.

trainIndex <- createDataPartition(dataset2$length, p=0.7, list=FALSE)
train <- dataset2[trainIndex,]
test <- dataset2[-trainIndex,]
testIndex <- createDataPartition(test$length, p=0.66, list=FALSE)
validate <- test[-testIndex,]
test <- test[testIndex,]

cat("train set: ", 100*nrow(train)/nrow(dataset2), "%", "\n", "test set: ", 100*nrow(test)/nrow(dataset2), "%", "\n", "validation set: ", 100*nrow(validate)/nrow(dataset2), "%")

## train set:  70.00494 % 
##  test set:  19.79955 % 
##  validation set:  10.1955 %

Regresja liniowa

Pierwszą wykorzystaną metodą będzie regresja liniowa. Do treningu zastosowano optymalizację według miary RMSE.

fit_lm <- train(length~., data=train, method="lm", metric="RMSE")

Obliczmy błąd popełniany przez model na zbiorze walidacyjnym w ujęciu RMSE oraz R^2.

predictions_lm <- predict(fit_lm, validate)
rmse_lm <- sqrt(mean((validate$length - predictions_lm)^2))
r2_lm <- cor(validate$length, predictions_lm) ^ 2
cat("RMSE: ", rmse_lm, "R2: ", r2_lm)

## RMSE:  1.362388 R2:  0.3292823

Regresja liniowa jest bardzo prostym modelem uczenia maszynowego o ograniczonych możliwościach, jednak błąd popełniony przez tą metodę nie jest bardzo duży.

SVM

Kolejnym regresorem wykorzystanym do predykcji długości śledzia jest metoda wektorów nośnych (ang. support vector machines, SVM). Implementacja tej metody znajduje się w pakiecie kernlab.

fit_svm <- ksvm(length~., train)
predictions_svm <- predict(fit_svm, validate)
rmse_svm <- sqrt(mean((validate$length - predictions_svm)^2))
r2_svm <- cor(validate$length, predictions_svm) ^ 2
cat("RMSE: ", rmse_svm, "R2: ", r2_svm)

## RMSE:  1.152424 R2:  0.5200524

Otrzymane wyniki są lepsze niż w przypadku regresji liniowej wedle miary RMSE, natomiast gorsze wedle R^2.W tym przypadku SVM trening trwał stosunkowo długo.

XGBoost

Ostatnią metodą jest eXtreme Gradient Boosting, w skrócie XGB. Wykorzystano implementację dostępną w pakiecie catboost. Metoda ta jest znana ze swojej wysokiej wydajności i efektywności. Wymaga innego przygotowania danych niż poprzednie metody - wymaga osobnego podania predyktorów i zmiennej wyjściowej, a także nie wymaga zakodowania one-hot zmiennych kategorialnych dzięki wbudowanemu mechanizmowi (konnieczność przekazania nazw zmiennych parametrem cat_features). Optymalizacja jest przeprowadzana domyslnie według miary RMSE.

trainIndex_xgb <- createDataPartition(dataset$length, p=0.7, list=FALSE)
train_xgb <- dataset[trainIndex_xgb,]
test_xgb <- dataset[-trainIndex_xgb,]
testIndex_xgb <- createDataPartition(test_xgb$length, p=0.66, list=FALSE)
validate_xgb <- test_xgb[-testIndex_xgb,]
test_xgb <- test_xgb[testIndex_xgb,]
train_Y <- train_xgb$length
test_Y <- test_xgb$length
validate_Y <- validate_xgb$length
train_xgb$length <- NULL
test_xgb$length <- NULL
validate_xgb$length <- NULL

train_pool <- catboost.load_pool(train_xgb, label=train_Y, cat_features = c("xmonth"))
test_pool <- catboost.load_pool(test_xgb, label=test_Y, cat_features = c("xmonth"))
val_pool <- catboost.load_pool(validate_xgb, label=validate_Y, cat_features = c("xmonth"))

fit_params <- list(iterations = 10000)
fit_catboost <- catboost.train(train_pool, test_pool = val_pool, params=fit_params)
predictions_xgb <- catboost.predict(fit_catboost, val_pool)
rmse_xgb <- sqrt(mean((validate$length - predictions_xgb)^2))
r2_xgb <- cor(validate$length, predictions_xgb) ^ 2

cat("RMSE: ", rmse_xgb, "R2: ", r2_xgb)

## RMSE:  1.591494 R2:  0.1752335

Uzyskane wyniki są najlepsze z dotychczas wykorzystanych metod uczenia maszynowego pod kątem miary R^2, jednak gorsze pod względem miary RMSE. Jednakże błąd R^2 jest prawie 2 razy mniejszy niż dla najlepszego dotychczasowego wyniku wedle tej miary, a RMSE tylko nieznacznie gorsze, dlatego dalszej analizie poddano tę metodę.

Na etapie analizy korelacji wyłoniono parę lcop1 i chel1 jako zmienne charakteryzujące sie najwyższym współczynnikiem korelacj równym 0.96. Spróbujmy usunąć zmienną chel1 i zaobserwować wpływ tego zabiegu na trafność predykcji modelu.

fit_params <- list(iterations = 10000, ignored_features=c(2))
fit_catboost_ch1 <- catboost.train(train_pool, test_pool = val_pool, params=fit_params)
predictions_xgb_ch1 <- catboost.predict(fit_catboost_ch1, val_pool)
rmse_xgb_ch1 <- sqrt(mean((validate$length - predictions_xgb_ch1)^2))
r2_xgb_ch1 <- cor(validate$length, predictions_xgb_ch1) ^ 2

cat("RMSE: ", rmse_xgb_ch1, "R2: ", r2_xgb_ch1)

## RMSE:  1.588176 R2:  0.1766239

Usunięcie chel2 ze zbioru atrybutów nie przyniosło znaczącej różnicy. Oznacza to jednak, że dla tej metody można wykorzystać mniej zmiennych bez utraty jakości predykcji, co jest dobrą informacją.

Analiza jakości predykcji na zbiorze testowym

Na koniec przetestowano zbudowane modele na zbiorze testowym.

test_lm <- predict(fit_lm, test)
test_svm <- predict(fit_svm, test)
test_xgb <- catboost.predict(fit_catboost, test_pool)
test_xgb_ch1 <- catboost.predict(fit_catboost_ch1, test_pool)

rmse_test_lm <- sqrt(mean((test$length - test_lm)^2))
r2_test_lm <- cor(test$length, test_lm) ^ 2
rmse_test_svm <- sqrt(mean((test$length - test_svm)^2))
r2_test_svm <- cor(test$length, test_svm) ^ 2
rmse_test_xgb <- sqrt(mean((test$length - test_xgb)^2))
r2_test_xgb <- cor(test$length, test_xgb) ^ 2
rmse_test_xgb_ch1 <- sqrt(mean((test$length - test_xgb_ch1)^2))
r2_test_xgb_ch1 <- cor(test$length, test_xgb_ch1) ^ 2

test_res <- data.frame("metoda" = c("regresja liniowa","SVM","XGB", "XGB bez chel1"), 
                "RMSE" = c(rmse_test_lm, rmse_test_svm, rmse_test_xgb, rmse_test_xgb_ch1), 
                "R2" =  c(r2_test_lm, r2_test_svm, r2_test_xgb, r2_test_xgb_ch1))

test_res %>%
  kable() %>%
  kable_styling()

Wyniki na zbiorze testowym wykazują podobne względne różnice w jakości predykcji jak dla zbioru walidacyjnego.

Analiza ważności atrybutów

Implementacja catboost pozwala na analizę ważności cech atrybutów. Przedstawia listę atrybutów wraz z ich procentową kontrybucją do ostatecznego wyniku predykcji.

catboost.get_feature_importance(fit_catboost)

##     cfin1     cfin2     chel1     chel2     lcop1     lcop2      fbar 
##  7.711926  7.286647  4.745681  5.103176  4.424734  5.706697  8.235664 
##      recr      cumf    totaln       sst       sal    xmonth       nao 
##  6.253652  5.489162  6.940089 17.395016  5.807484  6.963496  7.936575

Można zaobserwować, że najbardziej istotnym atrybutem okazała się temperatura przy powierzchni wody, która stanowi około 15% wkładu do predykcji.